Effects of Strong Correlations on the Zero Bias Anomaly in the Extended Hubbard 

Model with Disorder 
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We study the effect of strong correlations on the zero bias anomaly (ZBA) in disordered interacting 
systems. We focus on the two-dimensional extended Anderson-Hubbard model, which has both on- 
site and nearest-neighbor interactions on a square lattice. We use a variation of dynamical mean 
field theory in which the diagonal self-energy is solved self-consistently at each site on the lattice for 
each realization of the randomly-distributed disorder potential. Since the ZBA occurs in systems 
with both strong disorder and strong interactions, we use a simplified atomic-limit approximation 
for the diagonal inelastic self-energy that becomes exact in the large- disorder limit. The off-diagonal 
self-energy is treated within the Hartree-Fock approximation. The validity of these approximations 
is discussed in detail. We find that strong correlations have a significant effect on the ZBA at half 
filling, and enhance the Coulomb gap when the interaction is finite-ranged. 
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I. INTRODUCTION 

The term "zero bias anomaly" (ZBA) refers either to 
a peak in, or a suppression of, the density of states 
(DOS) at the Fermi energy. In disordered materials, a 
ZBA arises from the interplay between disorder and in- 
teractions. Zero bias anomalies were originally predicted 
to occur in strongly-disordered insulators by Efros and 
Shklovskii 1,2 (ES) and later in weakly-disordered met- 
als by Altshuler and Aronov 3 ( AA) . In the limit of weak 
interactions and disorder, AA showed that the exchange 
self-energy of the screened Coulomb interaction produces 
a cusp-like minimum in the DOS at the Fermi energy. In 
the limit of strong disorder, ES showed that the classical 
Hartree self-energy of the unscreened Coulomb interac- 
tion causes the DOS to vanish at the Fermi energy. Ex- 
periments have shown a smooth evolution between the 
AA and ES limits as a function of disorder^ and it ap- 
pears that the essential physics of the ZBA in conven- 
tional metals and insulators is well understood. 

In this work, we are interested in anomalies which have 
been observed in a number of transition metal oxide ma- 
terials, where the physics is less well understood. Tran- 
sition metal oxides often exhibit unconventional behav- 
ior because the physics of their valence band is domi- 
nated by strong short-ranged interactions whose effects 
cannot generally be explained by conventional theories 
of metals and insulators. Most notably, many transi- 
tion metal oxides exhibit a Mott transition when their 
valence band is half-filled. In disorder-free systems, the 
Mott transition^ occurs between a gapless metallic state 
and a gapped insulating state, and is driven by a strong 
intraorbital Coulomb interaction. 

The Mott transition may occur as a function of any 
number of parameters^ such as temperature^ or mag- 
netic field, 8 but more commonly occurs in transition- 
metal oxides as a result of chemical doping. A number 



of experiment o 9 ' 10 ! 11 ! 12 ? 13 have found that chemical dop- 
ing introduces sufficient disorder that there is a regime 
between the Mott-insulating and gapless phases that is 
characterized by a ZBA. This naturally raises the ques- 
tion of how the strong electron-electron correlations that 
are prevalent in the gapless phase near the Mott transi- 
tion affect the physics of the ZBA. 

The effect of strong correlations on the ZBA has re- 
ceived little attention^ In large part, this is due to the 
difficulty of incorporating both strong correlations and 
disorder in a manageable theory. A number of calcula- 
tions based on the unrestricted Hartree-Fock approxima- 
tion (HFA) have been used to study the phase diagram of 
the disordered Hubbard modeL 16 ' 17 i 18 ' 19 In these calcula- 
tions, the disorder potential is treated exactly for finite- 
sized systems, but the intraorbital Coulomb interaction 
is treated at the mean-field level and therefore neglects 
strong correlations. Much of the recent progress has in- 
volved various formulations of dynamical mean field the- 
ory (DMFT) to include disorder at some level of approx- 
imation. Coherent-potential-like approximations have 
been employed by a variety of authorsi 20 i 21 ' 22 i 23 ' 24 i 25 In 
these calculations, the local electron self-energy contains 
both inelastic contributions from the interactions and 
elastic contributions from the disorder-averaging process. 
It is well known that these kinds of disorder-averaging 
approximations capture many features of the DOS but 
do not retain the nonlocal correlations responsible for 
the ZBA and cannot, therefore, explain the experiments 
cited above^ An extension of DMFT, called statistical 
DMFT, has been employed to study ensembles consist- 
ing of Bethe lattices with random site energies! 26 ' 27 This 
represents an improvement over the disorder-averaged 
approximations in that the results depend nontrivially 
on the coordination number of the Bethe lattice. Very 
recently, the DMFT equations have been solved by us 
on a two dimensional square lattice in a way which pre- 
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serves spatial correlations between sitesi^ 8 . The calcula- 
tions employed a simple atomic-limit approximation for 
the self-energy that, while generally appropriate for the 
large-disorder limit, does not contain the off-diagonal 
self-energies responsible for the DOS anomalies observed 
in experiments. 

In this work, we extend our earlier calculations to in- 
clude a nearest neighbour interaction V. This interaction 
is treated at the mean-field level, and serves two pur- 
poses in this work: First, it represents finite range inter- 
actions which are present in real materials. Second, the 
exchange self-energy coming from this interaction plays 
a qualitatively similar role to the off-diagonal self-energy 
that is missing in our treatment of the intraorbital in- 
teraction U. The intraorbital self-energy is calculated 
within a Hubbard-I (HI) approximation which has the ef- 
fect of suppressing double occupancy of orbitals. Details 
of these calculations are outlined in Sec. Ill Al Because 
the HI approximation is known to be defficient in the 
disorder- free limit, we include a discussion in Sec. Ill Bl 
showing that the approximation for the diagonal self- 
energy is valid in the limit of large disorder and low co- 
ordination number. In Sec. Ill Ci we introduce a coherent 
potential approximation (CPA) for the disordered Hub- 
bard model which is used as a point of comparison for 
the lattice DMFT calculations. The results of numerical 
and analytical calculations are presented in Sec. IIIII 

Our primary result is that the ZBA is strongly en- 
hanced by electron correlations near half-filling when the 
interaction is finite-ranged. The ZBA is predominantly of 
the ES type, in the sense that the largest contribution is 
from the classical Hartree interaction between localized 
charges on neighboring sites. The enhancement of the 
ZBA at half-filling is due to strong correlations which in- 
hibit screening of the impurity potential: Near the Mott 
transition, the absence of screening drives the system to- 
wards the strongly localized limit where the physics of 
the Coulomb gap is most important. The magnitude of 
the ZBA is doping dependent, and drops off away from 
half-filling. 



II. CALCULATIONS 
A. Method 

We study an extended Anderson-Hubbard model, also 
known as a disordered t-U-V model. In this model, 
t denotes the electron kinetic energy, U the intraor- 
bital Coulomb interaction, and V the interaction be- 
tween nearest-neighbor sites. The V-term in the Hamil- 
tonian can be used to represent two distinct physical pro- 
cesses. First, it represents the nonlocal Coulomb interac- 
tion that, while neglected in the Hubbard model, is gen- 
erally present in real materials. Second, as shown below, 
the exchange self-energy from the nonlocal interaction 
is qualitatively similar to the exchange self-energy that 
arises in a perturbative treatment of the Hubbard model. 



This effective interaction plays a crucial role in clean low- 
dimensional systemsi 29 ' 30 ! 31 ' 32 It is often treated explic- 
itly in the large-£7 limit via an approximate mapping 
of the Hubbard model onto the t-J model, where J is 
the strength of the effective nonlocal interaction. A ma- 
jor difference between the t-J and extended Anderson- 
Hubbard models is that double occupation of orbitals is 
completely suppressed in the former whereas a finite frac- 
tion of sites will be doubly occupied in the latter when 
the width of the disorder distribution is larger than U. 
The extended Anderson-Hubbard model is 

H = -t c L c j<? + ^ 

(*,i> 

+ ^ ( e i^i + Un^hj{) , (1) 

i 

where denotes nearest-neighbor lattice sites i and 

j, hiu — ctj.Cjcr, hi — + hii, and parameters t, U 
and V are the kinetic energy, the on-site Coulomb inter- 
action, and the nearest-neighbor interaction respectively. 
ti is the site energy, which is box-distributed according 
to P(a) = W~ 1 Q{W/2 - \ei\), where W is the width of 
the disorder distribution and Q(x) the step-function. 

We treat the nearest-neighbor interaction at the mean- 
field level: 

/ \ ^ t 9 TijT1/q \ 

—n,n 3 w V \rnrij - ^ c i* c i°JH + fij Y~ J ^ 

with fji — (cLcif) = (CjiCn) in the paramagnetic phase 
and rij = (hj). Both fij and rij are determined self- 
consistently. The mean-field Hamiltonian, up to an ad- 
ditive constant, is: 

ft = + e ^ + U Y fli ^ flji ( 3 ) 

(i,j)a i i 

where t'^ = — t — V fij and e- = + VJ2j n ji where the 
sum is over nearest neighbors of i. 

The approximate Hamiltonian ^ is then solved using 
an iterative lattice DMFT (LDMFT) method that cap- 
tures the strong-correlation physics of the intraorbital 
interaction^ On an A^-site lattice, the single-particle 
Green's function can be expressed as an N X N matrix 
in the site-index: 

G(w) = [wl -t - e - EM]" 1 (4) 

with I the identity matrix, t the matrix of renormalized 
hopping amplitudes t'^ , e the diagonal matrix of renor- 
malized site energies e[ and the matrix of local 
self-energies. The self-energy £j(w) corresponds to the 
inelastic self-energy £(o>) (the so-called "impurity self- 
energy") in standard DMFT, which is obtained by vari- 
ous self-consistent impurity solvers* 3 ^ The iteration cycle 
begins with the calculation of G(oj) from Eq. ([5]), and 
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(5) 
(6) 



from the previous iteration. For each site i, one defines 
a Weiss mean field Gi(u>) = [G^cj) -1 + S^w)] -1 where 
Gij{uj) are the matrix elements of G(ui). The HI ap- 
proximation is the simplest improvement over the HFA 
that generates both upper and lower Hubbard bands 
and, as we discuss in the next section, it works well in 
the large-disorder limit W 3> t. In this approximation, 
= [G^uj)- 1 - Ef »] _1 where 



Y HI (< A - 77 ™ l + I 1 ' f ) 



(7) 



and rii is self-consistently determined for each site^ 

We remark that we can also express the Green's func- 
tion as 



G(w) = M-to-eo-SH]- 1 



(8) 



where to and eo are matrices of the unrenormalized hop- 
ping amplitudeds and site energies. In this case, the self- 
energy is: 



U- 



2 2 



and 



-vu 



(9a) 



(9b) 



This form emphasizes the nonlocal nature of the self- 
energy. 

B. Validity of the Self-Energy 

In this section, we discuss our treatment of the in- 
traorbital interaction in the renormalized Hamiltonian 
^ . We show that Eq. ((7|) is a reasonable approximation 
for the local self-energy in the large-disorder limit pro- 
vided that nonlocal effective interactions generated by U 
can be absorbed into the interaction V. 

Our discussion is based on a two-site Anderson- 
Hubbard Hamiltonian, ie. on Eq. ^ with two sites la- 
belled "1" and "2" . This Hamiltonian can, of course, 
be diagonalized exactly with relatively little effort. Here, 
we are interested in developing an approximate treatment 
that is valid in the large disorder limit, and which can be 
applied to the A-site problem. Comparison to the exact 
solution is used as a benchmark for the approximation. 



We use an equation-of-motion method to arrive at 
approximate expressions for the single-particle Green's 
function. Defining a Liouvillian superoperator C such 
that££ 



CA=[H,A], 



(10) 



where A is an arbitrary operator, we can formally write 
the time evolution of A as A(t) = cxp(iCt)A(0). It fol- 
lows directly that the retarded Green's function can be 
written 



G icr , ja (uj) = (41- 



1 



C 



A ) 



(11) 



where the inner product of two operators is defined as 
(A\B) = ({A 1 * , B}) and {, } refers to the anticommutator. 
The operator set c\ a is not closed under operations by 
£, but a complete operator set can be generated with 
repeated operation by £ on c\ a . For example, 

cci = 44, + 4Ai + u o>l - (is) 



where b\ a = ct^n^ — n^) and ~o = —a. The operator 
b\ a is a composite operator, and further composite oper- 
ators can be generated from C 2 c\ a , etc. The higher order 
composite operators involve excitations on multiple lat- 
tice sites, and are therefore expected to be less important 
in the disordered case than in the clean limit. Here, we 
truncate the series after a single application of C so that 
our operator basis consists of two operators, c il7 and b\ a , 
for each site and spin. This leads to a "two-pole" ap- 
proximation for the Green's function. This approach has 
been studied at length in the clean limit and has been 
shown to provide a reasonable qualitative description of 
the Hubbard model. 37 ! 38 ' 39 As shown in Fig. [TJ the two- 
pole approximation (described in more detail below) is 
essentially indistinguishable from the exact solution for 
the DOS of the two-site system. 

It is useful to define a generalized Green's function in 
the expanded operator space: 



M = 



f c t |_JL_ c t ) r c t I l fot ) 

f&t I 1 c t ) (bt l-^—b^- ) 



(13) 



such that Gi tT j a (u>) is given by the upper left quadrant 
of Giaj a (ui). Defining the Liouvillian matrix, 



(cUcbt 



(bl\Cc] a ) (bl\Cbl) 
and the matrix of overlap integrals, 



(14) 



/Cia.ja 



(4l<i) (4Jbl) 
(4l4) 0>3>l) 



(15) 
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FIG. 1: (Color online) Comparison of the approximate and 
exact densities of states for a two-site system. The site en- 
ergies are ei = W/2, 62 = — W/2 for W = St. Figure panels 
are for (a) U = it, (b) U = 8t and (c) U = 12t. Exact cal- 
culations are based on exact diagonalization of the two-site 
Hubbard model. Approximate densities of states are given by 
p(uj) — —\m^2 i Gii{u)/-K, with G(u>) the matrix (|19[) . Two 
different approximations are made for the self energies: the 
first has the self energy given by Eq. (|20p and makes the 
COMl approximation for p<22- the second has the self energy 
given by Eq. © with different choices for V. Both the COMl 
and exact solutions are offset for clarity. The COMl solution 
is essentially indistinguishable from the exact solution in all 
cases. Note that sf is the zero of energy in this and all other 
figures. 



we get 



with 



Gia,j a {u) = y/XivMiu - Hia\aVXj<y,ja (16) 



1 

ma(l - n m ) 



(17) 



and L = yj \ 1 L\/x 1 - 

For the nonmagnetic case, up and down spins 
are equivalent, and only the former need be consid- 
ered. For the two-site system, and taking the basis 
[cj-p cj-p b\f, b^], we write the Liouvillian matrix explic- 
itly as 



z\ + Un n -t' U x 

-t' e' 2 + Un 2l U 2 

Ui ei -t'p 

U 2 -t'p l 2 



(18) 



where 

Ui -- 

h - 

X - 

A. t = 

P z 
P : 



U^nnil-nn) 

e' q + U(l - nn) - t'Ai 
A,[n a (l - n^)}- 1 

((n iT - n^c^c^) - ((1 - n A - mi)A\ c i\) 

p[n n (l - n u )n 2L (l - n 2l )}~ 1/2 

(filial) - n n n 2l - (c T 1T c 2T (ct iCu + c\ L c 2l )) 



In the expression for A,, j is the nearest neighbor to i 
(i.e. j = 1 if i = 2 and j = 2 if i = 1). 

Equation [TH] allows us to solve for the Green's func- 
tion, provided the fields Aj, p, and mi are known. In 
practice, Aj and mi can be solved self-consistently, but 
further information is needed to calculate p. 

The Green's function G(u) is given by the upper left 
2x2 quadrant of Q(u>). It is straightforward to solve 
Eq. (fPo|) to show that 



G{w) = 



LO — 6 
t' 



! - En(w) 
- E 2 i(w) 



a - 

UJ — e' n 



E i2 (w) 
- E 22 (o>) 



with 



El2 



Unu + 



(w-ei)(w-e a )-t«pa 



-t'U 2 p 



(w-eO^-e^-t'V 



(19) 



(20a) 



(20b) 



where j is again the nearest neighbor to i. Equations 
|20a] and |20bj are the basic expressions for the self- 
energy. These expressions are shown in Fig. [T] to give 
very accurate results for the Green's function provided 
that p is correctly chosen. In this work, we have used 
the COMl approximation of Avella and Mancini. 39 The 
COMl approximation works well for the simple inhomo- 
geneous systems we have tested it on, but is extremely 
difficult to apply to disordered systems where p is differ- 
ent along every bond in the lattice. 

Some simplifications can be made in the case of large 
disorder. To begin with, we discuss the diagonal self- 
energies Ejj, and take i = 1. We note that p ~ 0(1), so 
that 



En ->■ Unii + 



U 2 n n (l-n n ) 



C/(l-n u )+i'Ai 



(21) 



whenever (u> — ei)(u> — e 2 ) ^S> t 2 . Apart from the shift 
tAi, this is just the Hubbard-I approximation for the 
self-energy. As we show next, Eq. [21] is justified for 
lu » Sf in the large disorder limit, which corresponds in 
the two-site problem to \ei — €2] ^S> t. 

We are interested in the validity of Eq. (f2T|) near ef 
and take the particular case £p = U/2 (which corre- 
sponds to half-filling) where strong correlations are most 
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important. When t = 0, each atomic Green's function 
has poles at e- and e' ; + U, so that spectral weight at Ep 
comes from sites with = ±E//2. This remains approx- 
imately true when t + 1 provided that W ^S> t. Then, if 
site 1 contributes spectral weight at the Fermi level, 



e F - ex 



U, 



and Eq. (|2"Tj) follows from Eq. (|20ap provided — e 2 | > 
t 2 /U. This condition will only not be met when e' 2 « 
— U/2. In other words, the two cases not well described 
by Eq. jUJ are (i) e' x » e£ and (ii) e' x « e' 2 + U. 

Certainly, in any random distributed set of site ener- 
gies, both cases are expected to occur for some fraction 
of sites on the lattice. However, if the disorder poten- 
tial is large, and the coordination number of the lattice 
is low, then the probability of any given site having a 
nearest neighbor satisfying either condition (i) or (ii) is 
low, and the fraction of sites not described by Eq. (f2~Tj) is 
small. It is interesting to note that the physical processes 
neglected here are (i) formation of singlet correlations be- 
tween nearly degenerate sites and (ii) resonant exchange 
between sites in which the double occupancy of site 1 is 
nearly degenerate with singlet formation between sites 1 
and 2. 

Two further comments are warranted regarding our 
treatment of En. First, the simplifications made above 
assume that both U and W are large. This case is di- 
rectly relevant to the current work since it is the regime in 
which density of states anomalies are observed. However, 
we have found empirically through numerical studies of 
small clusters that Eq. (|2ip is also a good approximation 
when U is small. This point is illustrated in Fig. QJa). 
Second, the term t'Ai in Eq. (f2~Tj) is neglected in Eq. (f9aj) . 
In the clean limit, |A| < 0.2 and depends only weakly 
on This term is small relative to the Hartree shift 
V J2j n j an d is therefore neglected. 

We emphasize that the Hartree shift V J2j n j does not 
arise naturally in our treatment of the Hubbard-[/ inter- 
action. It must therefore be undertstood to come directly 
from the nonlocal Coulomb interaction. This point is im- 
portant because, as is shown in Sec. lIIIl the Hartree term 
makes a large contribution to the ZBA at half-filling. 

The off-diagonal self-energy is significantly harder to 
evaluate than the diagonal term. It requires knowledge 
of a quantity p that cannot be calculated self-consistently 
within the two-pole approximation, although various ap- 
proximate schemes exist for its evaluation^ However, 
we note from the definition of p that it measures ex- 
change correlations between sites 1 and 2 and plays a 
similar role to the exchange self-energy — V/12 defined in 
Sec. Ill Al This suggests that, qualitatively, the exchange 
self-energy may represent the exchange term from the 
nonlocal Coulomb interaction. 

Figure[T]shows the DOS of a two-site system calculated 
within the approximation Q. The main effect of V is 
to produce a level repulsion between states above and 
below the Fermi energy, as illustrated in Fig. [IJb). The 
magnitude of the level repulsion decreases with increasing 



U, and is unobservably small for U — 12£ and V < St. 
As the figure shows, the DOS is well-reproduced with 
V = when U <C W and U ^> W, but is much less well 
reproduced when U ~ W. As expected from the analysis 
above, this can be corrected to some extent by a judicious 
choice of V. 



C. Coherent Potential Approximation 

In this section, we describe an implementation of the 
coherent potential approximation (CPA) that includes 
the effects of interactions and disorder within an effective 
medium approximation. As mentioned in the introduc- 
tion, the CPA neglects spatial correlations and therefore 
misses the physics of the ZBA. It is therefore useful as a 
point of comparison for our LDMFT calculations. 

Our CPA implementation applies specifically to the 
HFA and the HI approximation, and hence is denoted 
by HFA+CPA or HI+CPA as appropriate, and reduces 
to these approximations in the limit W — > 0. The 
HFA+CPA and HI+CPA algorithms also reduce to the 
usual CPA in the noninteracting U — > limit. Because 
of the local nature of these approximations, the nonlocal 
interaction V has not been included. 

For a particular lattice, we can calculate the local 
Green's function of the disorder-averaged system: 



GlocM 



N ^ 



1 



— €k — E(cj) 



(22) 



In this equation, is a self-energy that includes both 
inelastic contributions from the local interaction and elas- 
tic contributions from the disorder scattering. On the 
first iteration of the algorithm, E(w) is guessed, and on 
later iterations it is taken from the output of the previous 
iterations. Next, we define a Green's function 

G e (w) = [[G^Mr 1 + E(w) - e - E £ (u;)] ~* (23) 

This is the local Green's function for a site with energy 
e which is embedded in the effective medium. The term 
E e (<x>) is the inelastic self-energy for the site, and must 
be determined self-consistently. For both the HFA and 
HI approximation, E e (a>) depends on the local charge 
density n e . Equation ([23]) can therefore be closed by the 
relations 



2 f £F 
-— / Im G e (u>)dui, 
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and E e (w) = U^f for the HFA or 



2 uj-e-U(l-?f) 



(24) 



(25) 



for the HI approximation. Equations (j23j) -([25l) must be 
iterated to convergence for each value of e. 




FIG. 2: (Color online) Evolution of the DOS with doping for 
a 12 x 12 lattice with U = St, W = 12t, V = and 1000- 
2000 impurity configurations, (a)-(d) Black solid lines and 
red dashed lines represent the results of LDMFT and effective 
medium calculations respectively. The Lorentz broadening is 
7 = 0.025 throughout this work, (e)-(h) Corresponding plots 
for paramagnetic HFA calculations are shown for the same 
parameters and 1000 impurity configurations. 

We then average G e (u>) over site energies to get 

1 fW/2 

G av ( W ) = — / <feG e (w), (26) 

W J-W/2 

and a new self-energy is found via 

= [GtooH]" 1 + E(w) - [Gav^)]- 1 (27) 

The iteration cycle is now restarted at Eq. I|22p with 
E new (cj) taking the place of S(w). The iteration pro- 
cess is terminated when the difference between E new (u;) 
and £(w) is small. When the iteration cycle is complete, 
Gioc(w) is the disorder-averaged Green's function of the 
interacting system. 

III. RESULTS 

A. Numerical Results, V = 

We begin our discussion with the case V = 0. While it 
is not clear that a system can be experimentally realized 
in which the off-diagonal self-energy vanishes, this case 
is interesting because it provides a relatively simple illus- 
tration of the role of strong correlations. Throughout this 
section, strong correlation effects are identified by com- 
parisons between the HFA (which neglects correlations) 
and LDMFT. 

The DOS is calculated from the self-consistently deter- 
mined Green's function via 
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FIG. 3: Evolution of the DOS with U at quarter-filling for 
W = 12t, V = 0, and 8x8 lattices with 1000 impurity con- 
figurations. The upper panels (a)-(d) show LDMFT results 
while the lower panels (e)-(h) show HFA results. 



where the rank-iV matrix G(w) is given by Eq. The 
evolution of p(uS) with doping is shown for V = in 
Fig. O We have chosen W = I2t, which corresponds 
to W = 1.5D where D = 8t is the bandwidth in the 
clean noninteracting limit. We have taken U — 8t, which 
is large relative to t, but is much less than the critical 
U c ~ W at which the Mott transition takes place. For 
comparison, we have shown results for LDMFT and for 
the HI+CPA described in Sec. Iff CI These two theories 
employ the same approximation for the interaction and 
differ only in how they treat disorder: nonlocal spatial 
correlations between impurity sites and the charge den- 
sity are neglected in effective medium approximations 
and are treated exactly in LDMFT. The two methods 
give quantitatively similar results, except for the ZBA 
that emerges away from half filling in LDMFT but is 
absent in HI+CPA. 

We note that the sign of the ZBA in Fig. [5] is posi- 
tive. This is different from the usual case discussed in 
the literature, but is consistent with A A theory, where a 
negative ZBA comes from the exchange self-energy while 
the Hartree self-energy makes a weak positive correction 
to the DOS. Since the Hubbard interaction has a van- 
ishing exchange self-energy, the expectation from mean- 
field theory is for a positive ZBA when V = 0. This is 
illustrated by the numerical HFA calculations shown in 
Fig.lHe)-(h). 

Although the sign of the ZBA is the same in LDMFT 
and HFA calculations, its magnitude is different. In par- 
ticular, the peak at ep is finite at all doping levels in the 
HFA but is absent near half-filling in the LDMFT cal- 
culations. This difference shows that strong correlations 
suppress the ZBA near half filling. Results for quarter 
filling are shown in Fig. [3] for a range of U, where it can 
be seen that the ZBA grows with U when U is small, 
but saturates when U > 8t. A more technical discus- 



sion of these results is given in Sec. IHI CI and we briefly 
summarize the main ideas of this discussion here. 

The main distinction between weakly and strongly cor- 
related systems is that the local charge density m is a 
continuous variable in weakly correlated systems, but is 
restricted to near-integer values in strongly correlated 
systems. In the HFA, the energy of an isolated site is 
UJi = ei + Urii/2. For sites with Ep — U < u <£f, the self- 
consistent equation for the charge density, = 2f(u>i) 
[where f(x) is the Fermi function], will be satisfied at 
zero temperature by 

u> t = ef, Ui = 2(e F - £i)/U, 

where the second equality comes from rearranging the 
expression for u>i. Since a macroscopic fraction of sites 
satisfy sp ~U < e, < Ef, a peak (i.e. a positive ZBA) is 
expected in the DOS at the Fermi energy in the atomic 
limit. Numerical calculations (not shown) find that the 
peak persists, but weakens, as t/W grows. 

In contrast, the poles in the spectral function for an 
isolated strongly correlated site are at uii — ti and 

= £i + U. Because of the rigidity of the relationship be- 
tween UJi and e,, the distribution of uii values follows the 
distribution of and a vanishing fraction of sites there- 
fore have resonances at £p. There is, consequently, no 
ZBA when t = 0; the ZBA in LDMFT calculations only 
occurs when t/W is nonzero. The discussion in Sec. lIIICl 
shows that the spectral weight in the ZBA is proportional 
to the hybridization function Ai(£p ) between sites with 
Ci w Ef or ti + U ~ Ef and the rest of the lattice. The 
absence of a ZBA at half-filling comes from the fact that 
these sites decouple from the lattice, i.e. Ai(£p) = 0, when 
ef = U/2. 

Our analytical calculations in Sec. IIII CI also suggest 
that Ai(£p) is a strong function of both £p and U pro- 
vided £ f lies in the "central plateau" [by which we mean 
the broad peak in the DOS arising from the overlap 
of upper and lower Hubbard bands; see, for example, 
Fig-H^a)]. Outside the central plateau, Ai(£p) is a weak 
function of both £p and U . This is qualitatively con- 
sistent with the numerical results in Fig. [3j which show 
that the peak height increases with U for U < 8t and 
saturates at larger U : In the limit t/W — > 0, it is easy to 
show that ef lies in the central plateau for U < W/2. 

B. Numerical Results, V / 

We now consider the case V ^ 0. As for V = 0, 
the influence of strong correlations is most visible at half 
filling. For the purposes of this discussion, the term "ES- 
like behavior" refers to a negative Hartree contribution 
to the DOS, and the term "AA-like behavior" refers to a 
Hartree contribution which is positive and an exchange 
contribution which is negative. 

The DOS at half filling is shown in Fig. Hfor U = 8t 
and increasing V. In both the LDMFT and HFA results 
ES-like and AA-like behavior is present. As discussed 
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FIG. 4: Evolution of the DOS with V at half filling for U = St 
and W = 12t. (a)-(d) LDMFT results for an 8 x 8 lattice 
with 1000 impurity configurations, (e)-(h) HFA results are for 
the same parameters and 10 x 10 lattices with 1000 impurity 
configurations. 



in |40j there appears to be a transition as a function of 
energy, with ES-like behavior farther from the Fermi en- 
ergy and AA-like behavior closer in. Following the argu- 
ments of ES,— the average distance between states near 
the Fermi energy is very large. In our case the interac- 
tion has a finite range, and hence the ES behavior breaks 
down very near the Fermi energy. That the gap in Fig. [4] 
(d) is mostly ES-like, except right near the origin, can be 
seen in Fig. [5] (a) where the full result and a result with 
the Hartree term V J2j n j m Eq. EH set to zero are com- 
pared. Note that the DOS does not satisfy the ES result 
p{ui) oc \u> — Ep\ because the interaction is short-range. 
The most striking difference between the LDMFT and 
HFA results is the more pronounced ES-like behavior in 
the strongly-correlated case. Strong correlations result in 
much less screening of the disorder than in the mean-field 
treatment, hence enhancing the disorder-driven ES-like 
behavior. That the AA-like behaviors in the LDMFT 
and HFA results differ in sign is not surprising given the 
very different treatment of U in the two cases. 

The results at quarter filling are shown in Fig. [6] The 
ZBA in the LDMFT results is less pronounced at quar- 
ter filling than at half filling. The narrow ZBA close to 
the Fermi energy crosses over from positive to negative, 
consistent with increasing V and hence increasing nega- 
tive exchange contribution of the A A type. As seen in 
Fig. E](b), the Hartree self-energy modifies the DOS over 
a large energy range, but does not produce a gap-like fea- 
ture at the Fermi energy. Exact studies of small clusters, 
to be reported elsewhere^ suggest that the transition 
from large to small ZBA occurs when cf is shifted outside 
the central plateau described earlier. In the atomic limit, 
when e p lies within the central plateau rii may have three 
distinct values (0,1 or 2), whereas rtj can only be or 1 
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FIG. 5: (color online) Contributions to the zero bias anomaly 
for V = 2.4t at (a) half filling and (b) quarter filling. For 
the curve labeled "Full" , both the exchange and Hartree self- 
energy contributions from the nearest neighbor interaction 
are retained, while for the curve labelled "Exchange Only", 
the Hartree self-energy is set to zero. The curves show that 
the ZBA comes primarily from the Hartree contribution, and 
is therefore of the Efros-Shklovskii type. The inset shows 
results for the localization length at half-filling reproduced 
from a finite-size scaling analysis in Ref. [2S| for V = 0. The 
short localization length at large U in LDMFT calculations, 
relative to HF calculations, is consistent with the enhanced 
Coulomb gap in the LDMFT calculations. 



for tp below the lower edge of the plateau. In small clus- 
ters, this reduction in the range of possible charge states 
for individual sites directly results in a reduced ZBA. In 
marked contrast to the strongly-correlated results, the 
HFA results show an AA-like peak and an ES-like dip, 
both of which grow with increasing V. 

Figure [7] shows the variation of the DOS with U for 
V = 1.6t. For the somewhat artificial case of U = 
shown in Fig. 6(a) and (e), the results differ slightly be- 
cause the numerical HFA and LDMFT routines converge 
differently. Both obtain charge-density-wave order frus- 
trated by the disorder, but with slightly shifted domain 
walls. In the HFA case, increasing U screens the disor- 
der, reducing the ES-like behavior. The AA-like Hartree 
peak is initially strengthened by U but is reduced as the 
screening increases. At large U [Fig. [7Jh)], the DOS ap- 
praches the clean-limit result. In the LDMFT case, the 
screening produced by U initially weakens the ES-like 
behavior. However, for larger values of U, the screen- 
ing in the strongly correlated case is much less than that 
in mean field. The localization length based on a finite- 
size scaling analysis (calculated for V = 0) is reproduced 
from a previous work2& in the inset of Fig. [5] While the 
localization length grows monotonically with U in the 
HFA, in the LDMFT it reaches a maximum at U ~ At 
and decreases with increasing U thereafter. Similarly, the 
ES-like behavior of the LDMFT DOS is initially weak- 
ened, but is not lost as in the HFA and saturates before 
the opening of the Mott gap. 

We note that, in Fig. the exchange contribution to 
the ZBA depends only weakly on doping. We showed pre- 
viously in Sec. Ill Bl that the exchange self-energy arising 
from the nonlocal interaction plays the same role as the 
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FIG. 6: Evolution of the DOS with V at quarter filling with 
W = 12t and U = St. Calculations are for an 8 x 8 site 
lattice with more than 1000 sample configurations for each 
parameter set. Results are shown for (a)-(d) LDFMT and 
(e)-(h) HFA. 
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FIG. 7: Evolution of the DOS with U at half-filling with 
W = Ylt and V = l.Gt. (a)-(d) LDMFT and (e)-(h) HFA 
calculations are shown. Calculations are for an 8 x 8 site 
lattice with more than 1000 sample configurations for each 
parameter set. 



exchange self-energy arising from a perturbative treat- 
ment of intraorbital interaction. Thus, the exchange-only 
curves shown in Fig. [5] should behave the same qualita- 
tively as an exact treatment of the V = Anderson- 
Hubbard model. Such a study has been reported by 
Chiesa et al.^ where the ZBA was found to be indepen- 
dent of doping for fillings 0.6 < n < 1. This is consistent 
with our findings. 

When the Hartree self-energy is included, our findings 
differ substantially from Chiesa et al. In particular, the 
doping-depcndence of the ZBA due to ES-like physics 
which we have reported here is not found when the inter- 
action is purely locals We emphasize that the Hartree 
self-energy responsible for the ZBA is physically distinct 
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from the exchange contributions arising from the Hub- 
bard U interaction, and can only occur when there is 
a finite-ranged interaction. In summary, the difference 
in results between Chiesa et al. and us appears to stem 
from the difference in the range of the electron-electron 
interactions. 



C. Analysis of strong correlation effects on the 
zero bias anomaly 

In this section, we discuss the origin of the ZBA in 
the V = case. We begin with Eq. ([4]) for the Green's 
function. In the large disorder limit, it is possible to treat 
the hopping matrix element as a perturbation. When the 
matrix t is zero, G(uj) decouples into a diagonal matrix 
describing an ensemble of isolated atoms with Green's 
functions 



1 



(29) 



where the superscript zeros refer to the isolated atomic 
systems and E° (ui) is given exactly by Eq. ((?]). The 
atomic Green's function, G° (w), has poles at 



with spectral weights 



ZL = l - 



u. 



7 o _ "» 
i+ ~ ~2 



(30) 



(31) 



The total density of states is found by averaging the 
imaginary part of G% (u>) over e, in the interval — W/2 < 
ti < W/2. Since the pole energies are linear functions of 
ej, the total density of states is featureless at the Fermi 
energy. 

Next, we use the fact that the diagonal matrix elements 
of Eq. ^ can be written in the form 



1 



■ Aj(w) - £i(w) 



(32) 



where Ai(u>) is the hybridization function that describes 
the coupling of site i to the rest of the lattice. It is 



(33) 



where Uj are the matrix elements of t between sites i 

and j, G^(lo) is the Green's function for the lattice with 
site i removed, and the second line is the expansion of the 
first to 0(t 2 ). Equation ([33"]) applies formally to the limit 
that the localization length vanishes, but is qualitatively 
correct for t -C W. We note that Aj(w) is a complex 
function of frequency for metallic systems, but is real 
with a discrete spectrum of simple poles for Anderson- 
localized systems as we have here. 



Recalling Eq. |[7J), we solve for the poles of Gu{uS) to 

o(t 2 y. 



Ai(w;_), 



(34a) 
(34b) 



where we have assumed Aj <C U . In the approximation 
(|33p . Ai(uji±) diverges when Ui± is degenerate with any 
u>j± for nearest neighbor site j. This is an artifact of the 
approximation since any degeneracy between i and j is 
lifted by hybridization of the orbitals. The poles of Ai(u) 
must therefore differ from ui%± by an energy > t, and we 
impose a cutoff |Aj(cj$±)| < t. 

The spectral weights Zi± of the poles (|3"4"|) are reduced 
by 0(t 2 ) from Zf± and the remaining spectral weight 
appears at new poles resulting from hybridization of site 
i with the rest of the lattice. These poles play a role in 
suppressing the ZBA in the limit that the localization 
length becomes large, but are of secondary importance 
when t -C W as in the current discussion. 

Equations (|3"4"1) contain the essential physics of the 
ZBA, which we summarize here before we go into the 
detailed calculations. In both equations, the local charge 
susceptibility \a — ~drii/dei is nonzero because of the 
hybridization function Aj(w). The main idea is that, be- 
cause Xu is nonzero, sites with energies et that are suf- 
ficiently close to Ef (ejt — U) can adjust their filling rii 
such that LOi- = £f {^>i+ = £f)- The range of ti satisfy- 
ing the criterion of "sufficiently close" is set by Aj(£p), 
and the weight under the ZBA peak is therefore also set 
by Ai(ep). The suppression of the ZBA at half filling 
then follows from the fact that the disorder average of 
Ai(£f) is an antisymmetric function of Ef- 

We consider sites with energies ej such that u>i± = Ef- 
The expression for u>i± requires knowledge of the charge 
density n,, which is given by rii/2 = J2± ^i±f( UJ i±) + 
0{t 2 ). For sites with ua± — £f, this reduces to [excepting 
terms of 0(t 2 )] 



Hi. 

~2 



(l-y)/M, (n i <l) 



(35) 



for Ui- = ef and 



HI 

2 



l-y)+y/(£F), K>1) (36) 



for tOi+ — Ef- At zero temperature, < /(ef) < 1 and 
these equations are satisfied for a range of e*. Setting 
uji- = ef in Eq. (I34a|) and applying the restriction < 
Ui < 1, we generate the limits el < £f - £i < £[; on e„ 
where 



(-ii 



i„m,'^U( e ,; 



2 

Ai(ef) 



,Ai(e F ) 



In this range [rearranging Eq. (|34ap ] 

_ Ui _ E F ~ £» 

2 ~ Ai(E F ) ' 



(37) 
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and the density of states at sf coming from sites with 
u)i- = £p is therefore 

5p_(uKE F ) = — / deiZi^5(uj - e F ) 



3 |A(ef)1 
8 W 



5{uj-e F ) + 0(t i ). (38) 



In this equation, |A(£f)| is an average over |Aj(£f)|- An 
identical result can be found for sites with resonance en- 
ergies u>i + = sf, so that the total density of states near 
the Fermi energy is 



p(u ~ E F ) 



3 [AM] 

4 W 



S(uj — e F ). 



(39) 



Equation (|39[) shows that the ZBA is a delta- function at 
zero temperature. At finite temperatures T, the ZBA is 
a peak of width ~ T. The spectral weight in the peak 
is proportional to the hybridization function A(e F ), and 
we next show how this depends on doping. 

We consider a site i in a lattice with coordination 
number Z c , whose nearest neighbors have randomly cho- 
sen site energies. At half filling (e F = U/2), the terms 
GjJsF) in the sum in ([33| are positive or negative with 
equal probability, and tend to cancel. As the filling is re- 
duced, the probability that G^fep) is negative (positive) 
becomes larger (smaller). We make a rough calculation 
that illustrates this behavior by replacing the sum over 
site index j in (|3"3")) with an integral over Cj . Thus 



t. Equation (|40aP applies when the Fermi level sits in the 
central plateau, and shows that Ai(e F ) is antisymmetric 
about half- filling (i.e. e F = U/2), and grows linearly away 
from half- filling. Outside of the central plateau, Ai(e F ) 
is a weak function of e F . 

In order to compare with Fig. [31 we evaluate Eq. (|40b[) 
at quarter filling, which for small t corresponds to 



e F 



U _ W jj w 

2 ±> a, 

n TT \ w 



Then Eq. (|40b)) gives 



|A*M| 



t 2 z c 

W 




(W_y_ip 
W -U 



2 

















for U < \ , and 



\Hsf)\ 



t 2 Z c 
2W 



(41a) 
(41b) 



for U > -f . In Eqs. gTj, |Ai(e F )| grows linearly with 
U for small U (recall that there is a cutoff such that 
ln(U/t) — > \n(t/t) when U < t), and saturates at a finite 
value when U ^> Both these results, and the results 
at half filling in Eq. (|40aj) are qualitatively consistent 
with the numerical results shown in Figs. [2]and[3l 



A 4 (£p) 



t 2 Z c 

w 



W/2 



de 



W/2 



1 - n e /2 



n e /2 



e F 



e F — e — U 



where n e is the charge density for sites with site-energy 
e. Recalling our constraint |Ai(ei?)| < t, we introduce 
cutoffs near the poles of the integrand. Noting that 

2, e<e F -U 
1, e F — U < e < e F 
0, e > e F 



we get 



A l {e F ) 



t 2 z c 
w 



In 



(40a) 



for U - ¥f < e F < ¥ and 



A,(e F ) 



In 



^ + s F 
~w 

— -e F 



In 



U-%--e F 



w i _ 

— +£f 



■ih.fl 

2 V ti 



(40b) 



for 



< C7 



— . The logarithmic divergences in 



— < e F - - 2 
Eqs. (|40p are artificial and must be cut off whenever any 
numerator or denominator has a magnitude smaller than 



IV. CONCLUSIONS 

We have studied the effects of strong correlations on 
the zero bias anomaly in the density of states for disor- 
dered interacting systems. Our results show significant 
doping dependence. When only local interactions are in- 
cluded, a positive ZBA is suppressed by strong correla- 
tions at half filling due to cancelation in the hybridiza- 
tion function responsible for the peak. When nearest- 
neighbor interactions are included (simultaneously pro- 
ducing a more accurate treatment of the Hubbard U 
term), strong correlations modify the mean field results 
at both half and quarter filling. In particular, at half fill- 
ing ES-like behavior is enhanced due to reduced screening 
in the strongly correlated system, whereas at quarter fill- 
ing the ZBA is smaller reflecting the doping dependence 
of the ES physics. 
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